Lipid Analysis

R Markdown

load the libraries

library("readxl")
library("tidyverse")
library("rstatix")
library("DESeq2")
library("EnhancedVolcano")
library("ggplot2")
library("pmp")
library("ggpubr")
library("ComplexHeatmap")
library("limma")
library("qgraph")

Dataset-1

load the files and assign them to respective variable names.

lipid <-as.data.frame(read_excel("italian cohort lipid study.xlsx",sheet = 2)) 

lipid<-lipid[-1,-277]

clinical_data<-lipid[,1:4]

lipid_data<-lipid[,5:ncol(lipid)]

annotations <- as.data.frame(read_excel("annotations.xlsx",col_names =T))

colnames(lipid_data)<-make.unique(annotations$lipid)

rownames(lipid_data)<-clinical_data$sample_ID

lipid_data<-as.data.frame(lipid_data)

summary(lipid_data[,1:10])

##     CE_14:0             CE_15:0             CE_16:0            CE_16:0.1       
##  Min.   :0.0000000   Min.   :1.885e-05   Min.   :0.0001391   Min.   :0.003984  
##  1st Qu.:0.0003600   1st Qu.:8.394e-05   1st Qu.:0.0002549   1st Qu.:0.006530  
##  Median :0.0004367   Median :1.192e-04   Median :0.0004512   Median :0.007532  
##  Mean   :0.0004403   Mean   :1.187e-04   Mean   :0.0005237   Mean   :0.007500  
##  3rd Qu.:0.0005176   3rd Qu.:1.476e-04   3rd Qu.:0.0006468   3rd Qu.:0.008625  
##  Max.   :0.0008344   Max.   :2.630e-04   Max.   :0.0016558   Max.   :0.010846  
##     CE_16:1            CE_16:1.1           CE_17:0             CE_17:1         
##  Min.   :0.000e+00   Min.   :0.001046   Min.   :2.787e-05   Min.   :0.0000000  
##  1st Qu.:3.230e-05   1st Qu.:0.002081   1st Qu.:7.312e-05   1st Qu.:0.0001557  
##  Median :5.871e-05   Median :0.002471   Median :1.012e-04   Median :0.0002392  
##  Mean   :8.938e-05   Mean   :0.002615   Mean   :1.105e-04   Mean   :0.0002002  
##  3rd Qu.:1.196e-04   3rd Qu.:0.003015   3rd Qu.:1.453e-04   3rd Qu.:0.0002709  
##  Max.   :4.621e-04   Max.   :0.006832   Max.   :2.295e-04   Max.   :0.0003762  
##     CE_18:1            CE_18:1.1      
##  Min.   :0.0003304   Min.   :0.01170  
##  1st Qu.:0.0005593   1st Qu.:0.02253  
##  Median :0.0008584   Median :0.02588  
##  Mean   :0.0009193   Mean   :0.02585  
##  3rd Qu.:0.0011769   3rd Qu.:0.02856  
##  Max.   :0.0023625   Max.   :0.04234

table(clinical_data$NAFLD)

## 
##   0   1 
## 120  21

scale the data

lipid_scaled<-scale(lipid_data)

summary(lipid_scaled[,1:10])

##     CE_14:0            CE_15:0            CE_16:0          CE_16:0.1       
##  Min.   :-3.11746   Min.   :-2.21915   Min.   :-1.0770   Min.   :-2.35097  
##  1st Qu.:-0.56875   1st Qu.:-0.77224   1st Qu.:-0.7526   1st Qu.:-0.64815  
##  Median :-0.02568   Median : 0.01217   Median :-0.2030   Median : 0.02164  
##  Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.54740   3rd Qu.: 0.64303   3rd Qu.: 0.3445   3rd Qu.: 0.75223  
##  Max.   : 2.79007   Max.   : 3.20715   Max.   : 3.1696   Max.   : 2.23800  
##     CE_16:1          CE_16:1.1          CE_17:0           CE_17:1       
##  Min.   :-1.0552   Min.   :-1.9408   Min.   :-1.7820   Min.   :-1.8814  
##  1st Qu.:-0.6738   1st Qu.:-0.6612   1st Qu.:-0.8060   1st Qu.:-0.4180  
##  Median :-0.3621   Median :-0.1781   Median :-0.2014   Median : 0.3673  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.0000  
##  3rd Qu.: 0.3566   3rd Qu.: 0.4948   3rd Qu.: 0.7513   3rd Qu.: 0.6645  
##  Max.   : 4.3995   Max.   : 5.2148   Max.   : 2.5666   Max.   : 1.6547  
##     CE_18:1          CE_18:1.1        
##  Min.   :-1.3858   Min.   :-2.662815  
##  1st Qu.:-0.8471   1st Qu.:-0.625705  
##  Median :-0.1432   Median : 0.004171  
##  Mean   : 0.0000   Mean   : 0.000000  
##  3rd Qu.: 0.6064   3rd Qu.: 0.509735  
##  Max.   : 3.3961   Max.   : 3.101539

conduct PCA analysis

res_pca_1<-prcomp(lipid_scaled,centre=T,scale=F)

res_pca <- as.data.frame(res_pca_1$x)

res_pca$steatosis<-clinical_data$NAFLD

res_pca$sample<-clinical_data$sample_ID


theme<-theme(panel.background = element_blank(),panel.border=element_rect(fill=NA),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),strip.background=element_blank(),axis.text.x=element_text(colour="black"),axis.text.y=element_text(colour="black"),axis.ticks=element_line(colour="black"),plot.margin=unit(c(1,1,1,1),"line"))


percentage <- round(res_pca_1$sdev / sum(res_pca_1$sdev) * 100, 2)
percentage <- paste( colnames(res_pca), "(", paste( as.character(percentage), "%", ")", sep="") )


plot1<-ggplot(res_pca,aes(x=PC1,y=PC2,label= sample,color=steatosis ))

plot1<-plot1+geom_point(size=4)+theme+ xlab(percentage[1]) + ylab(percentage[2])

print(plot1)

rename the sample in clinical data

clinical_data$NAFLD[clinical_data$NAFLD==0]<-"steatosis0"

clinical_data$NAFLD[clinical_data$NAFLD==1]<-"steatosis1"

head(clinical_data)

##   sample_ID patient_no      NAFLD sex
## 2        N1          1 steatosis0   F
## 3        N2          1 steatosis0   F
## 4        N3          1 steatosis0   F
## 5        N4          1 steatosis0   F
## 6        N5          1 steatosis0   F
## 7        N6          1 steatosis0   F

Extract the steatosis 0 and steatosis 1 sample separately

clinical_data_steatosis0<-subset(clinical_data, NAFLD=="steatosis0",select=sample_ID:sex)

clinical_data_steatosis1<-subset(clinical_data, NAFLD=="steatosis1",select=sample_ID:sex)

head(clinical_data_steatosis1)

##    sample_ID patient_no      NAFLD sex
## 76       N75         15 steatosis1   F
## 77       N76         15 steatosis1   F
## 82       N81         18 steatosis1   F
## 83       N82         18 steatosis1   F
## 86       N85         20 steatosis1   M
## 87       N86         20 steatosis1   M

head(clinical_data_steatosis0)

##   sample_ID patient_no      NAFLD sex
## 2        N1          1 steatosis0   F
## 3        N2          1 steatosis0   F
## 4        N3          1 steatosis0   F
## 5        N4          1 steatosis0   F
## 6        N5          1 steatosis0   F
## 7        N6          1 steatosis0   F

lets see how many steatosis0 and steatosis1 samples are present

nrow(clinical_data_steatosis0)

## [1] 120

nrow(clinical_data_steatosis1)

## [1] 21

There are more steatosis0 than steatosis 1 , hence lets conduct differential expression of lipid analysis between these two groups by subseting the samples into separate batches.

create 6 batches of distinct (without sample replacement of steatosis0) and conduct differential expression of lipid for each batch of steatosis0 and steatosis1 such as (20 vs 21) and identify the repeatedly arising significant gene between them.

Lets create 6 batches

rownames(clinical_data_steatosis0)<-1:120

batch1_clinical_data_steatosis0<-clinical_data_steatosis0[1:20,]
  
batch2_clinical_data_steatosis0<-clinical_data_steatosis0[21:40,]

batch3_clinical_data_steatosis0<-clinical_data_steatosis0[41:60,]

batch4_clinical_data_steatosis0<-clinical_data_steatosis0[61:80,]

batch5_clinical_data_steatosis0<-clinical_data_steatosis0[81:100,]

batch6_clinical_data_steatosis0<-clinical_data_steatosis0[101:120,]

now create the clinical data for those 6 batches by merging clinical batch of steatosis0 and steatosis1

batch1_clinical_data<-rbind(batch1_clinical_data_steatosis0,clinical_data_steatosis1)

batch2_clinical_data<-rbind(batch2_clinical_data_steatosis0,clinical_data_steatosis1)

batch3_clinical_data<-rbind(batch3_clinical_data_steatosis0,clinical_data_steatosis1)

batch4_clinical_data<-rbind(batch4_clinical_data_steatosis0,clinical_data_steatosis1)

batch5_clinical_data<-rbind(batch5_clinical_data_steatosis0,clinical_data_steatosis1)

batch6_clinical_data<-rbind(batch6_clinical_data_steatosis0,clinical_data_steatosis1)

extract the expression data for each of the batches now

batch1_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch1_clinical_data$sample_ID,]

batch2_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch2_clinical_data$sample_ID,]

batch3_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch3_clinical_data$sample_ID,]

batch4_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch4_clinical_data$sample_ID,]

batch5_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch5_clinical_data$sample_ID,]

batch6_lipid_scaled<-lipid_scaled[rownames(lipid_scaled) %in% batch6_clinical_data$sample_ID,]

Differential lipid analysis

conduct differential lipid expression analysis for batch-1

design_batch1 <- model.matrix(~0+batch1_clinical_data$NAFLD,ref="steatosis0")

colnames(design_batch1) <- c("Steatosis0","Steatosis1")

lipid_data_limma1<-t(batch1_lipid_scaled)

fit_data1 <- lmFit(lipid_data_limma1, design_batch1)

contrasts1 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch1)

fit2_data1 <- contrasts.fit(fit_data1, contrasts1)

fit2_data1<- eBayes(fit2_data1)

full_results1 <- topTable(fit2_data1, number=Inf,adjust = "BH")

full_results1 <- cbind(symbols = rownames(full_results1), full_results1)

rownames(full_results1) <- 1:nrow(full_results1)

lets view the plot for batch 1 and extract the significant genes of them

p_cutoff <- 0.05
fc_cutoff <- 1
topN <-10

full_results1 %>% 
  mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>% 
  mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>% 
  ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

significant_batch1<-full_results1[full_results1$adj.P.Val<p_cutoff,]

head(significant_batch1)

##   symbols     logFC      AveExpr         t      P.Value    adj.P.Val        B
## 1 TG_53:3  1.781418  0.007213676  10.84643 1.509766e-14 3.931815e-12 22.91345
## 2 DG_34:3  2.125923 -0.021706774  10.56533 3.714068e-14 3.931815e-12 22.01915
## 3 DG_43:6 -1.832345  0.029401360 -10.51652 4.346753e-14 3.931815e-12 21.86286
## 4 TG_56:4  2.004968  0.186773092  10.33753 7.759063e-14 3.931815e-12 21.28707
## 5 TG_56:5  2.036917  0.106174432  10.33529 7.815784e-14 3.931815e-12 21.27983
## 6 TG_50:3  1.766674 -0.181929281  10.30776 8.547425e-14 3.931815e-12 21.19090

dim(significant_batch1)

## [1] 191   7

conduct differential lipid expression analysis for batch-2

design_batch2 <- model.matrix(~0+batch2_clinical_data$NAFLD,ref="steatosis0")

colnames(design_batch2) <- c("Steatosis0","Steatosis1")

lipid_data_limma2<-t(batch2_lipid_scaled)

fit_data2 <- lmFit(lipid_data_limma2, design_batch2)

contrasts2 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch2)

fit2_data2 <- contrasts.fit(fit_data2, contrasts2)

fit2_data2<- eBayes(fit2_data2)

full_results2 <- topTable(fit2_data2, number=Inf,adjust = "BH")

full_results2 <- cbind(symbols = rownames(full_results2), full_results2)

rownames(full_results2) <- 1:nrow(full_results2)

volcano plot for batch 2

full_results2 %>% 
  mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>% 
  mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>% 
  ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

significant_batch2<-full_results2[full_results2$adj.P.Val<p_cutoff,]

head(significant_batch2)

##   symbols     logFC     AveExpr          t      P.Value    adj.P.Val        B
## 1 DG_38:5 -1.833957 0.008044289 -11.115877 3.355189e-15 9.260321e-13 24.37448
## 2 DG_43:6 -1.834214 0.030313308 -10.684875 1.383043e-14 1.908599e-12 22.96750
## 3 PI_39:4 -1.693211 0.032347451 -10.212007 6.717783e-14 6.180360e-12 21.39697
## 4 PS-O_38 -1.738166 0.021359607 -10.010747 1.327170e-13 9.157471e-12 20.72026
## 5 TG_56:4  1.928321 0.224161721   9.829341 2.461660e-13 1.358836e-11 20.10622
## 6 TG_53:3  1.766265 0.014605727   9.568558 6.022766e-13 2.770473e-11 19.21690

dim(significant_batch2)

## [1] 189   7

conduct differential lipid expression for batch 3

design_batch3 <- model.matrix(~0+batch3_clinical_data$NAFLD,ref="steatosis0")

colnames(design_batch3) <- c("Steatosis0","Steatosis1")

lipid_data_limma3<-t(batch3_lipid_scaled)

fit_data3 <- lmFit(lipid_data_limma3, design_batch3)

contrasts3 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch3)

fit2_data3 <- contrasts.fit(fit_data3, contrasts3)

fit2_data3<- eBayes(fit2_data3)

full_results3 <- topTable(fit2_data3, number=Inf,adjust = "BH")

full_results3 <- cbind(symbols = rownames(full_results3), full_results3)

rownames(full_results3) <- 1:nrow(full_results3)

volcano plot for batch 3

full_results3 %>% 
  mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>% 
  mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>% 
  ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

significant_batch3<-full_results2[full_results3$adj.P.Val<p_cutoff,]

head(significant_batch3)

##   symbols     logFC     AveExpr          t      P.Value    adj.P.Val        B
## 1 DG_38:5 -1.833957 0.008044289 -11.115877 3.355189e-15 9.260321e-13 24.37448
## 2 DG_43:6 -1.834214 0.030313308 -10.684875 1.383043e-14 1.908599e-12 22.96750
## 3 PI_39:4 -1.693211 0.032347451 -10.212007 6.717783e-14 6.180360e-12 21.39697
## 4 PS-O_38 -1.738166 0.021359607 -10.010747 1.327170e-13 9.157471e-12 20.72026
## 5 TG_56:4  1.928321 0.224161721   9.829341 2.461660e-13 1.358836e-11 20.10622
## 6 TG_53:3  1.766265 0.014605727   9.568558 6.022766e-13 2.770473e-11 19.21690

dim(significant_batch3)

## [1] 171   7

conduct differential lipid expression for batch 4.

design_batch4 <- model.matrix(~0+batch4_clinical_data$NAFLD,ref="steatosis0")

colnames(design_batch4) <- c("Steatosis0","Steatosis1")

lipid_data_limma4<-t(batch4_lipid_scaled)

fit_data4 <- lmFit(lipid_data_limma4, design_batch4)

contrasts4 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch4)

fit2_data4 <- contrasts.fit(fit_data4, contrasts4)

fit2_data4<- eBayes(fit2_data4)

full_results4 <- topTable(fit2_data4, number=Inf,adjust = "BH")

full_results4 <- cbind(symbols = rownames(full_results4), full_results4)

rownames(full_results4) <- 1:nrow(full_results4)

volcano plot for batch 4

full_results4 %>% 
  mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>% 
  mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>% 
  ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

significant_batch4<-full_results4[full_results4$adj.P.Val<p_cutoff,]

head(significant_batch4)

##   symbols     logFC     AveExpr         t      P.Value    adj.P.Val         B
## 1 DG_38:5 -1.583231 -0.22742195 -7.799571 4.108759e-10 5.608425e-08 12.923874
## 2 DG_34:3  1.536189  0.47264808  7.771243 4.539084e-10 5.608425e-08 12.826047
## 3 DG_43:6 -1.569167 -0.22818895 -7.687448 6.096114e-10 5.608425e-08 12.536394
## 4 PS-O_38 -1.269659 -0.29150250 -7.349281 2.012102e-09 1.388351e-07 11.363847
## 5 Oxochol  1.285460 -0.05608858  6.472756 4.516807e-08 2.493277e-06  8.312214
## 6 PI_41:6 -1.219338 -0.22827882 -6.095570 1.720613e-07 7.914820e-06  7.003015

dim(significant_batch4)

## [1] 129   7

Conduct differential lipid expression for batch 5

design_batch5 <- model.matrix(~0+batch5_clinical_data$NAFLD,ref="steatosis0")

colnames(design_batch5) <- c("Steatosis0","Steatosis1")

lipid_data_limma5<-t(batch5_lipid_scaled)

fit_data5 <- lmFit(lipid_data_limma5, design_batch5)

contrasts5 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch5)

fit2_data5 <- contrasts.fit(fit_data5, contrasts5)

fit2_data5<- eBayes(fit2_data5)

full_results5 <- topTable(fit2_data5, number=Inf,adjust = "BH")

full_results5 <- cbind(symbols = rownames(full_results5), full_results5)

rownames(full_results5) <- 1:nrow(full_results5)

volcano plot for batch 5

full_results5 %>% 
  mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>% 
  mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>% 
  ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

significant_batch5<-full_results5[full_results5$adj.P.Val<p_cutoff,]

head(significant_batch5)

##     symbols      logFC     AveExpr         t      P.Value    adj.P.Val        B
## 1 DG_36:2.1 -1.1993450  0.49101262 -5.958626 4.239642e-07 6.882647e-05 6.304671
## 2   PC_36:4  1.5397202  0.08691141  5.910031 4.987425e-07 6.882647e-05 6.147818
## 3   TG_53:2 -0.9889126  0.87228183 -4.635981 3.318158e-05 2.184353e-03 2.109665
## 4   PC_38:4  1.1610951 -0.12761206  4.623782 3.450870e-05 2.184353e-03 2.072151
## 5   TG_53:3 -0.7102970  0.86506896 -4.510687 4.957409e-05 2.184353e-03 1.725896
## 6   TG_53:6 -1.2651848  0.69742118 -4.465661 5.722572e-05 2.184353e-03 1.588845

dim(significant_batch5)

## [1] 42  7

Conduct different lipid expression for batch 6

design_batch6 <- model.matrix(~0+batch6_clinical_data$NAFLD,ref="steatosis0")

colnames(design_batch6) <- c("Steatosis0","Steatosis1")

lipid_data_limma6<-t(batch6_lipid_scaled)

fit_data6 <- lmFit(lipid_data_limma6, design_batch6)

contrasts6 <- makeContrasts(Steatosis1 - Steatosis0, levels=design_batch6)

fit2_data6 <- contrasts.fit(fit_data6, contrasts6)

fit2_data6 <- eBayes(fit2_data6)

full_results6 <- topTable(fit2_data6, number=Inf,adjust = "BH")

full_results6 <- cbind(symbols = rownames(full_results6), full_results6)

rownames(full_results6) <- 1:nrow(full_results6)

volcano plot for batch 6

full_results6 %>% 
  mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>% 
  mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>% 
  ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

significant_batch6<-full_results6[full_results6$adj.P.Val<p_cutoff,]

head(significant_batch6)

##   symbols      logFC     AveExpr         t      P.Value    adj.P.Val        B
## 1 PC-O_44  1.4429320  0.10057341  5.973130 4.210383e-07 0.0001069434 6.326011
## 2 TG_54:2 -0.9965229  0.21923572 -5.685326 1.095100e-06 0.0001069434 5.405205
## 3 PC-O_36  1.3393702 -0.06098412  5.667329 1.162428e-06 0.0001069434 5.347761
## 4 PS-P_28  0.9521793 -0.14752453  5.376447 3.041048e-06 0.0002098323 4.422556
## 5 PC-O_40  1.4239131  0.06407383  5.038354 9.218303e-06 0.0005088503 3.357679
## 6 DG_36:2 -1.0665274  0.37025704 -4.825684 1.839361e-05 0.0008461062 2.695810

dim(significant_batch6)

## [1] 57  7

Find the intersecting significant lipids among the 6 batches

lipids_batch1<-significant_batch1$symbols

lipids_batch2<-significant_batch2$symbols

lipids_batch3<-significant_batch3$symbols

lipids_batch4<-significant_batch4$symbols

lipids_batch5<-significant_batch5$symbols

lipids_batch6<-significant_batch6$symbols



intersected_lipids<-Reduce(intersect,list(lipids_batch1,lipids_batch2,lipids_batch3,lipids_batch4,lipids_batch5,lipids_batch6))

print(intersected_lipids)

##  [1] "TG_53:3"     "DG_34:1"     "TG_52:3"     "TG_52:2.1"   "TG_53:2"    
##  [6] "TG_53:7"     "TG_54:2"     "DG_36:2"     "TG_53:6"     "CE_16:0.1"  
## [11] "Cholesterol"

conduct different expression of lipid only with the intersected lipids

common_lipid_scaled<-lipid_scaled[,colnames(lipid_scaled) %in% intersected_lipids]

design_final <- model.matrix(~0+clinical_data$NAFLD,ref="steatosis0")

colnames(design_final) <- c("Steatosis1","Steatosis0")

lipid_data_limma_final<-t(common_lipid_scaled)

fit_data_final <- lmFit(lipid_data_limma_final, design_final)

contrasts_final <- makeContrasts(Steatosis1 - Steatosis0, levels=design_final)

fit2_data_final <- contrasts.fit(fit_data_final, contrasts_final)

fit2_data_final<- eBayes(fit2_data_final)

full_results_final <- topTable(fit2_data_final, number=Inf,adjust = "BH")

full_results_final <- cbind(symbols = rownames(full_results_final), full_results_final)

rownames(full_results_final) <- 1:nrow(full_results_final)

print(full_results_final)

##        symbols      logFC       AveExpr         t      P.Value    adj.P.Val
## 1      TG_53:3 -1.0295330  3.430077e-17 -4.531666 6.307155e-06 0.0000693787
## 2    CE_16:0.1  0.9155465  2.225613e-16  4.029935 5.853538e-05 0.0003219446
## 3      TG_52:3 -0.8818559 -5.731354e-17 -3.881640 1.081696e-04 0.0003904006
## 4      TG_53:7 -0.8665732  5.344425e-17 -3.814371 1.419638e-04 0.0003904006
## 5      DG_36:2 -0.8525358  4.763723e-17 -3.752583 1.815657e-04 0.0003994444
## 6    TG_52:2.1 -0.8391009  7.714598e-17 -3.693447 2.290206e-04 0.0004198711
## 7      TG_53:2 -0.7702227  3.874216e-17 -3.390267 7.160239e-04 0.0011251804
## 8      DG_34:1 -0.7479340 -1.896877e-16 -3.292160 1.016891e-03 0.0013982253
## 9      TG_54:2 -0.7253484 -4.478293e-17 -3.192745 1.437954e-03 0.0017574991
## 10 Cholesterol  0.6866663  3.254021e-16  3.022480 2.548805e-03 0.0028036859
## 11     TG_53:6 -0.4913213  1.675669e-17 -2.162635 3.072386e-02 0.0307238648
##             B
## 1   3.5407333
## 2   1.4405388
## 3   0.8669213
## 4   0.6138069
## 5   0.3852174
## 6   0.1699332
## 7  -0.8800780
## 8  -1.2006081
## 9  -1.5158089
## 10 -2.0331984
## 11 -4.2129728

visualize

full_results_final %>% 
  mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff ) %>% 
  mutate(Rank = 1:n(), Label = ifelse(Rank < topN,symbols,"")) %>% 
  ggplot(aes(x = logFC, y = -log10(P.Value), col=Significant,label=Label)) + geom_point() + geom_text_repel(col="black")+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

significant_final<-full_results_final[full_results_final$adj.P.Val<p_cutoff,]

head(significant_final)

##     symbols      logFC       AveExpr         t      P.Value    adj.P.Val
## 1   TG_53:3 -1.0295330  3.430077e-17 -4.531666 6.307155e-06 0.0000693787
## 2 CE_16:0.1  0.9155465  2.225613e-16  4.029935 5.853538e-05 0.0003219446
## 3   TG_52:3 -0.8818559 -5.731354e-17 -3.881640 1.081696e-04 0.0003904006
## 4   TG_53:7 -0.8665732  5.344425e-17 -3.814371 1.419638e-04 0.0003904006
## 5   DG_36:2 -0.8525358  4.763723e-17 -3.752583 1.815657e-04 0.0003994444
## 6 TG_52:2.1 -0.8391009  7.714598e-17 -3.693447 2.290206e-04 0.0004198711
##           B
## 1 3.5407333
## 2 1.4405388
## 3 0.8669213
## 4 0.6138069
## 5 0.3852174
## 6 0.1699332

dim(significant_final)

## [1] 11  7

lets observe the -log10(P.value) of the significant lipids with a plot

graph<-ggplot(data=significant_final, aes(x=symbols, y=-log10(P.Value))) +
  geom_bar(stat="identity", fill="steelblue")+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))

graph

observe the lipids logFC of the lipid

graph_FC<-ggplot(data=significant_final, aes(x=symbols, y=logFC)) +
  geom_bar(stat="identity", fill="steelblue")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"))

graph_FC

Non-paramteric test

perform wilcoxon test on the significant lipids

#extract the significant lipids

lipid_stat<-lipid_scaled[,colnames(lipid_scaled) %in% significant_final$symbols ]

number.of.samples<-ncol(lipid_stat)

# Create a output

p.value.vector <- c()

# Assign it to a list

p.value.list <- numeric(number.of.samples)

# Add the sample type and preprocess the dataframe

lipid_stat<-as.data.frame(cbind(lipid_stat,group=clinical_data$NAFLD))

lipid_steatosis0<-lipid_stat[lipid_stat$group=="steatosis0", ]

lipid_steatosis1<-lipid_stat[lipid_stat$group=="steatosis1", ]

lipid_steatosis0<-t(lipid_steatosis0)

lipid_steatosis1<-t(lipid_steatosis1)

input.combined.b <- cbind(lipid_steatosis0[,1:120],lipid_steatosis1[,1:21])

input.combined.b<-as.data.frame(input.combined.b)

input.combined.b<-input.combined.b[-12,]

input.combined.b = data.frame(lapply(input.combined.b, function(x) as.numeric(as.character(x))),check.names=F, row.names = rownames(input.combined.b))


# Establishing variable with number of observations
number.of.samples <-nrow(input.combined.b)


# Establishing vector with number of observations
vector.number.of.samples <-1:number.of.samples

for (i in vector.number.of.samples) {
  steatosis0 <- unlist(input.combined.b[i, 1:120])
  steatosis1 <- unlist(input.combined.b[i, 121:141])
  wilcox.test.i <- wilcox.test(steatosis0, steatosis1,
                               mu=0,
                               alt="two.sided",
                               paired=F,
                               conf.int=F,
                               conf.level=0.95,
                               exact=F
  )
  p.value.list[i] <- wilcox.test.i["p.value"]
}

create a dataframe

#lipid list names to make a summary

lipid.list.df <- paste(rownames(input.combined.b) )

lipid.list.df<-as.data.frame(lipid.list.df)
p.value<- data.frame(matrix(unlist(p.value.list ), nrow=length(p.value.list ), byrow=TRUE))

colnames(p.value)<-"p.value"

#binding lipid list and pvalue to create summary list

wilcox.test.results <- cbind(lipid.list.df, p.value)

print(wilcox.test.results)

##    lipid.list.df      p.value
## 1      CE_16:0.1 7.555428e-05
## 2    Cholesterol 1.376790e-03
## 3        DG_34:1 5.172790e-04
## 4        DG_36:2 9.612430e-05
## 5      TG_52:2.1 5.364438e-05
## 6        TG_52:3 3.013631e-05
## 7        TG_53:2 1.296171e-04
## 8        TG_53:3 2.821378e-06
## 9        TG_53:6 3.473401e-03
## 10       TG_53:7 2.849056e-05
## 11       TG_54:2 8.942429e-05

print(wilcox.test.results$lipid.list.df)

##  [1] "CE_16:0.1"   "Cholesterol" "DG_34:1"     "DG_36:2"     "TG_52:2.1"  
##  [6] "TG_52:3"     "TG_53:2"     "TG_53:3"     "TG_53:6"     "TG_53:7"    
## [11] "TG_54:2"

visualize the p values from wilxocon test

bar_graph<-ggplot(data=wilcox.test.results, aes(x=lipid.list.df, y=-log10(p.value))) +
  geom_bar(stat="identity", fill="steelblue")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
 
#theme_minimal()

bar_graph

Order the clinical data

clinical_data_boxplot<-clinical_data[match(colnames(input.combined.b), clinical_data$sample_ID),]

boxplot_data<-t(input.combined.b)

boxplot_data<-as.data.frame(boxplot_data)

boxplot_data<-add_column(boxplot_data, samples = clinical_data_boxplot$NAFLD, .before = 1)

head(boxplot_data)

##       samples  CE_16:0.1 Cholesterol   DG_34:1    DG_36:2  TG_52:2.1   TG_52:3
## N1 steatosis0 1.16329901  0.18285982 -1.251192 -1.2132063 -0.8970818 -1.222416
## N2 steatosis0 0.43946354 -0.19793041 -1.246836 -0.8620915 -0.8583501 -1.315292
## N3 steatosis0 0.03958326 -0.30009357 -1.452855 -1.2032769 -0.9457813 -1.262850
## N4 steatosis0 0.16038153 -0.05679797 -1.336211 -0.8620326 -0.9107814 -1.199429
## N5 steatosis0 0.33983494  0.13088503 -1.318653 -0.8693971 -0.8836531 -1.241908
## N6 steatosis0 0.49945675  0.73867030 -1.389051 -1.3103400 -0.9648729 -1.208322
##       TG_53:2    TG_53:3   TG_53:6   TG_53:7    TG_54:2
## N1 -1.2597033 -1.2412499 -1.251368 -1.057828 -0.9565674
## N2 -0.9259357 -0.5980203 -1.251368 -1.057828 -0.7283276
## N3 -1.2597033 -0.9276351 -1.251368 -1.057828 -0.9942757
## N4 -1.0089002 -1.2412499 -1.251368 -1.057828 -0.7656261
## N5 -0.9301686 -0.8299763 -1.251368 -1.004041 -0.7705726
## N6 -1.2597033 -0.9003478 -1.251368 -1.057828 -0.8283395

#write.csv(boxplot_data,"lipids_for_analysis.csv")

create a boxplot for visualizing the difference between the steatosis-0 and steatosis-1 samples.

boxplot_data$samples=as.factor(boxplot_data$samples)

my_lipid_list <- c("CE_16:0.1","Cholesterol","DG_34:1","DG_36:2","TG_52:2.1","TG_52:3","TG_53:2","TG_53:3","TG_53:6","TG_53:7","TG_54:2")

my_plot_list <- vector(mode = "list", length = 11)

for (i in 1:11) {
  my_comparisons <- list( c("steatosis0", "steatosis1") )
  p <-ggboxplot(boxplot_data, x = "samples", y = my_lipid_list[i],
          color = "samples", palette = "jco",add = "jitter")+
          stat_compare_means(label.y = 16)
  my_plot_list[[i]] <- p
}
 
print(my_plot_list[1:11])

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

Visualization - Heatmap

create complex heatmap with all the samples in a single plot

ht1 = Heatmap(input.combined.b[,1:120], name = "steatosis0", row_title = "", column_title = "steatosis0")
ht2 = Heatmap(input.combined.b[121:141], name = "steatosis1", row_title = "", column_title = "steatosis1")
ht_list = ht1 + ht2

draw(ht_list, row_title = "Lipids", row_title_gp = gpar(col = "black"),
     column_title = "Comparison between steatosis0 and steatosis1", column_title_side = "bottom")

lets annotate the samples along with their gender to view their differences in heatmap

clinical_data_F<-subset(clinical_data, sex=="F",select=sample_ID:sex)

clinical_data_M<-subset(clinical_data, sex=="M",select=sample_ID:sex)

#lets add their gender 

clinical_data_F$sample_sex<-sub("^","F.",clinical_data_F$sample_ID)

clinical_data_M$sample_sex<-sub("^","M.",clinical_data_M$sample_ID)

clinical_data_both_sex<-rbind(clinical_data_F,clinical_data_M)

#order the newly created metadata to the colnames of the lipid expression data

clinical_data_both_sex_ordered<-clinical_data_both_sex[match(colnames(input.combined.b), clinical_data_both_sex$sample_ID),]

lets check if the column names of the expression lipid data is similar to the clinical data

all(colnames(input.combined.b) == clinical_data_both_sex_ordered$sample_ID)

## [1] TRUE

head(input.combined.b[,1:10])

##                     N1         N2          N3          N4         N5         N6
## CE_16:0.1    1.1632990  0.4394635  0.03958326  0.16038153  0.3398349  0.4994568
## Cholesterol  0.1828598 -0.1979304 -0.30009357 -0.05679797  0.1308850  0.7386703
## DG_34:1     -1.2511917 -1.2468364 -1.45285479 -1.33621146 -1.3186531 -1.3890511
## DG_36:2     -1.2132063 -0.8620915 -1.20327690 -0.86203257 -0.8693971 -1.3103400
## TG_52:2.1   -0.8970818 -0.8583501 -0.94578129 -0.91078143 -0.8836531 -0.9648729
## TG_52:3     -1.2224155 -1.3152917 -1.26284959 -1.19942874 -1.2419081 -1.2083217
##                     N7         N8          N9         N10
## CE_16:0.1   -0.7594497 -1.0764997 -0.50037377  0.06703728
## Cholesterol -0.9287608 -0.8771608 -0.33348685  0.08801765
## DG_34:1     -1.1093491 -0.5819656 -0.80144507 -1.26994037
## DG_36:2     -1.1479058  0.3503971 -0.79397421 -0.97467791
## TG_52:2.1   -0.8771772 -0.6298855 -0.75281917 -0.92615116
## TG_52:3     -0.9137006 -0.7715755 -0.08312642 -1.22462884

now lets replace the column names of the expression data with the newly created sample_ID+gender column

colnames(input.combined.b)<-clinical_data_both_sex_ordered$sample_sex

head(input.combined.b[,1:10])

##                   F.N1       F.N2        F.N3        F.N4       F.N5       F.N6
## CE_16:0.1    1.1632990  0.4394635  0.03958326  0.16038153  0.3398349  0.4994568
## Cholesterol  0.1828598 -0.1979304 -0.30009357 -0.05679797  0.1308850  0.7386703
## DG_34:1     -1.2511917 -1.2468364 -1.45285479 -1.33621146 -1.3186531 -1.3890511
## DG_36:2     -1.2132063 -0.8620915 -1.20327690 -0.86203257 -0.8693971 -1.3103400
## TG_52:2.1   -0.8970818 -0.8583501 -0.94578129 -0.91078143 -0.8836531 -0.9648729
## TG_52:3     -1.2224155 -1.3152917 -1.26284959 -1.19942874 -1.2419081 -1.2083217
##                   F.N7       F.N8        F.N9       F.N10
## CE_16:0.1   -0.7594497 -1.0764997 -0.50037377  0.06703728
## Cholesterol -0.9287608 -0.8771608 -0.33348685  0.08801765
## DG_34:1     -1.1093491 -0.5819656 -0.80144507 -1.26994037
## DG_36:2     -1.1479058  0.3503971 -0.79397421 -0.97467791
## TG_52:2.1   -0.8771772 -0.6298855 -0.75281917 -0.92615116
## TG_52:3     -0.9137006 -0.7715755 -0.08312642 -1.22462884

extract the file used to create the heatmap

file_heatmap<-t(input.combined.b)

file_heatmap<-as.data.frame(file_heatmap)

all(rownames(file_heatmap) == clinical_data_both_sex_ordered$sample_sex)

## [1] TRUE

file_heatmap<-add_column(file_heatmap,sample=clinical_data_both_sex_ordered$NAFLD,.before=1)

#head(file_heatmap)
#dim(file_heatmap)

#write.csv(file_heatmap,"lipid_analysis.csv")

view the heatmap for significant lipids by subsetting the samples into smaller categories.

ht1 = Heatmap(input.combined.b[,1:30], name = "steatosis0", row_title = "", column_title = "steatosis0")

ht2 = Heatmap(input.combined.b[,31:60], name = "steatosis0", row_title = "", column_title = "steatosis0")

ht3 = Heatmap(input.combined.b[,61:90], name = "steatosis0", row_title = "", column_title = "steatosis0")

ht4 = Heatmap(input.combined.b[,91:120], name = "steatosis0", row_title = "", column_title = "steatosis0")

draw(ht1, row_title = "Lipids", row_title_gp = gpar(col = "black"),
     column_title = " steatosis0", column_title_side = "bottom")

draw(ht2, row_title = "Lipids", row_title_gp = gpar(col = "black"),
     column_title = " steatosis0 ", column_title_side = "bottom")

draw(ht3, row_title = "Lipids", row_title_gp = gpar(col = "black"),
     column_title = " steatosis0", column_title_side = "bottom")

draw(ht4, row_title = "Lipids", row_title_gp = gpar(col = "black"),
     column_title = " steatosis0", column_title_side = "bottom")

ht5 = Heatmap(input.combined.b[,121:141], name = "steatosis1", row_title = "", column_title = "steatosis1")

draw(ht5, row_title = "Lipids", row_title_gp = gpar(col = "black"),
     column_title = " steatosis1", column_title_side = "bottom")

lets create a correlation network

set.seed(111)

Data2<-t(input.combined.b)

corMat <- cor(Data2) # Correlate data

head(corMat)

##              CE_16:0.1 Cholesterol    DG_34:1    DG_36:2  TG_52:2.1    TG_52:3
## CE_16:0.1    1.0000000   0.8500356 -0.7254860 -0.7240950 -0.7718605 -0.5241541
## Cholesterol  0.8500356   1.0000000 -0.5952215 -0.6494509 -0.6631440 -0.3506181
## DG_34:1     -0.7254860  -0.5952215  1.0000000  0.8576533  0.9495881  0.7917609
## DG_36:2     -0.7240950  -0.6494509  0.8576533  1.0000000  0.8757690  0.6876598
## TG_52:2.1   -0.7718605  -0.6631440  0.9495881  0.8757690  1.0000000  0.7471021
## TG_52:3     -0.5241541  -0.3506181  0.7917609  0.6876598  0.7471021  1.0000000
##                TG_53:2    TG_53:3    TG_53:6    TG_53:7    TG_54:2
## CE_16:0.1   -0.6755248 -0.7026092 -0.5763073 -0.6052691 -0.6623299
## Cholesterol -0.5533712 -0.5618185 -0.4802842 -0.4557635 -0.5734958
## DG_34:1      0.8932375  0.8427915  0.8708539  0.8297891  0.8073444
## DG_36:2      0.7601000  0.7698806  0.7926969  0.7282854  0.8003767
## TG_52:2.1    0.9030812  0.8679832  0.8465218  0.8285776  0.8019836
## TG_52:3      0.7377178  0.7564427  0.8649190  0.9224849  0.6230012

plot the graph

Groups <- c("Cholesteryl ester","Cholesterol",rep("Diacylglycerol",2),rep("Triacylglycerol",7))

Graph_pcor <- qgraph(corMat, graph = "pcor", layout = "spring",theme = "colorblind", threshold = "BH",
                     sampleSize = nrow(Data2), alpha = 0.05,groups=Groups,
                     legend.cex = 0.55, vsize = 10,
                    esize = 15, borders = T,label.fill.vertical=0.3,
                    label.scale=T,edge.labels=T,edge.label.bg=T,repulsion=0.35)

lets obtain the CSV files of the significant data for the random forest model

#write.csv(clinical_data,'clinical_data_significant.csv')

#write.csv(input.combined.b,'lipid_data_significant.csv')

obtain the session info

sessionInfo()

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
##  [1] grid      stats4    parallel  stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] igraph_1.2.6                qgraph_1.6.9               
##  [3] limma_3.46.0                ComplexHeatmap_2.6.2       
##  [5] ggpubr_0.4.0                pmp_1.2.1                  
##  [7] EnhancedVolcano_1.8.0       ggrepel_0.9.1              
##  [9] DESeq2_1.30.1               SummarizedExperiment_1.20.0
## [11] Biobase_2.50.0              MatrixGenerics_1.2.1       
## [13] matrixStats_0.59.0          GenomicRanges_1.42.0       
## [15] GenomeInfoDb_1.26.7         IRanges_2.24.1             
## [17] S4Vectors_0.28.1            BiocGenerics_0.36.1        
## [19] rstatix_0.7.0               forcats_0.5.1              
## [21] stringr_1.4.0               dplyr_1.0.6                
## [23] purrr_0.3.4                 readr_1.4.0                
## [25] tidyr_1.1.3                 tibble_3.1.1               
## [27] ggplot2_3.3.3               tidyverse_1.3.1            
## [29] readxl_1.3.1               
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.1             tidyselect_1.1.1       RSQLite_2.2.7         
##   [4] AnnotationDbi_1.52.0   htmlwidgets_1.5.3      BiocParallel_1.24.1   
##   [7] munsell_0.5.0          codetools_0.2-16       withr_2.4.2           
##  [10] colorspace_2.0-1       highr_0.9              ggalt_0.4.0           
##  [13] knitr_1.33             rstudioapi_0.13        ggsignif_0.6.1        
##  [16] Rttf2pt1_1.3.8         GenomeInfoDbData_1.2.4 mnormt_2.0.2          
##  [19] bit64_4.0.5            vctrs_0.3.8            generics_0.1.0        
##  [22] xfun_0.23              itertools_0.1-3        randomForest_4.6-14   
##  [25] R6_2.5.0               ggbeeswarm_0.6.0       clue_0.3-59           
##  [28] locfit_1.5-9.4         bitops_1.0-7           cachem_1.0.5          
##  [31] DelayedArray_0.16.3    assertthat_0.2.1       scales_1.1.1          
##  [34] nnet_7.3-14            beeswarm_0.4.0         gtable_0.3.0          
##  [37] ash_1.0-15             Cairo_1.5-12.2         rlang_0.4.10          
##  [40] genefilter_1.72.1      GlobalOptions_0.1.2    splines_4.0.3         
##  [43] extrafontdb_1.0        impute_1.64.0          broom_0.7.6           
##  [46] checkmate_2.0.0        yaml_2.2.1             reshape2_1.4.4        
##  [49] abind_1.4-5            modelr_0.1.8           backports_1.2.1       
##  [52] Hmisc_4.5-0            extrafont_0.17         tools_4.0.3           
##  [55] lavaan_0.6-9           psych_2.1.6            ellipsis_0.3.2        
##  [58] jquerylib_0.1.4        RColorBrewer_1.1-2     Rcpp_1.0.6            
##  [61] plyr_1.8.6             base64enc_0.1-3        zlibbioc_1.36.0       
##  [64] RCurl_1.98-1.3         rpart_4.1-15           pbapply_1.4-3         
##  [67] GetoptLong_1.0.5       haven_2.4.1            cluster_2.1.0         
##  [70] fs_1.5.0               magrittr_2.0.1         data.table_1.14.0     
##  [73] openxlsx_4.2.3         circlize_0.4.13        reprex_2.0.0          
##  [76] pcaMethods_1.82.0      tmvnsim_1.0-2          hms_1.1.0             
##  [79] evaluate_0.14          xtable_1.8-4           XML_3.99-0.6          
##  [82] rio_0.5.26             jpeg_0.1-8.1           gridExtra_2.3         
##  [85] shape_1.4.6            compiler_4.0.3         maps_3.3.0            
##  [88] KernSmooth_2.23-17     crayon_1.4.1           htmltools_0.5.1.1     
##  [91] corpcor_1.6.9          Formula_1.2-4          geneplotter_1.68.0    
##  [94] lubridate_1.7.10       DBI_1.1.1              dbplyr_2.1.1          
##  [97] proj4_1.0-10.1         MASS_7.3-53            Matrix_1.2-18         
## [100] car_3.0-10             cli_2.5.0              pkgconfig_2.0.3       
## [103] foreign_0.8-80         xml2_1.3.2             foreach_1.5.1         
## [106] pbivnorm_0.6.0         annotate_1.68.0        vipor_0.4.5           
## [109] bslib_0.2.5.1          missForest_1.4         XVector_0.30.0        
## [112] rvest_1.0.0            digest_0.6.27          rmarkdown_2.8         
## [115] cellranger_1.1.0       htmlTable_2.2.1        curl_4.3.1            
## [118] gtools_3.8.2           rjson_0.2.20           glasso_1.11           
## [121] lifecycle_1.0.0        nlme_3.1-149           jsonlite_1.7.2        
## [124] carData_3.0-4          fansi_0.4.2            pillar_1.6.1          
## [127] lattice_0.20-41        ggrastr_0.2.3          fastmap_1.1.0         
## [130] httr_1.4.2             survival_3.2-7         glue_1.4.2            
## [133] zip_2.1.1              fdrtool_1.2.16         png_0.1-7             
## [136] iterators_1.0.13       bit_4.0.4              stringi_1.5.3         
## [139] sass_0.4.0             blob_1.2.1             latticeExtra_0.6-29   
## [142] memoise_2.0.0